from loguru import logger
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
import sep
from astropy.wcs import WCS
from datetime import datetime, timedelta
# some useful functions

import numpy as np
from astropy.stats import sigma_clip, sigma_clipped_stats
from astropy.modeling import models, fitting
import matplotlib.pyplot as plt
import os, sys
from   collections import Counter
from   scipy.spatial import cKDTree as ckdt
import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord, search_around_sky
from collections import Counter
from scipy.spatial import KDTree
import os


def todegree(ra_hms, dec_dms):
    from astropy.coordinates import Angle
    import astropy.units as u

    # Convert RA (in hours, minutes, seconds) to degrees
    ra_angle = Angle(ra_hms, unit=u.hourangle)
    ra_deg = ra_angle.deg

    # Convert Dec (in degrees, minutes, seconds) to degrees
    dec_angle = Angle(dec_dms, unit=u.deg)
    dec_deg = dec_angle.deg
    return [ra_deg, dec_deg]


def calculate_limiting_magnitude(image_path, catalog_path, ra_dec, band, mag_zp, target_region_half_size, output_path, filter_,PIXSCAL):
    """
    计算给定图像的极限星等 (3σ 和 5σ)，并可视化目标区域。

    参数：
        image_path (str): 输入图像文件路径，FITS 格式。
        catalog_path (str): 星表文件路径（在当前版本中似乎未实际使用）。
        ra_dec (tuple): 目标天体的 (RA, Dec) 坐标，用于找到对应的像素位置。
        band (str): 滤光片名称，比如 'u','v','g','r','i','z','y' 等。
        mag_zp (float): 零点，用于将测得的亮度转换成星等。
        target_region_half_size (float): 目标区域的一半大小(单位为角秒)。
        output_path (str): 输出图像（可视化）保存路径。
        filter_ (str): 滤光片名称，用于命名输出文件。
    
    返回：
        limiting_mag_info (list): 包含图像名称、滤光片、曝光时间、FWHM、背景 RMS、3σ/5σ 极限星等、以及零点的列表。
    """
    # 读取图像数据（FITS 格式），并提取数据、头部信息、曝光时间以及 WCS 坐标信息
    with fits.open(image_path) as hdul:
        img_data_non_scaled = hdul[0].data
        header = hdul[0].header
        exp_time = header['EXPTIME']  # 获取曝光时间
        wcs = WCS(header)             # 从头部信息初始化 WCS

    # 将 RA, Dec 转换为图像中的像素坐标（浮点），再取整
    pixel_positions = wcs.world_to_pixel_values(ra_dec[0], ra_dec[1])
    pixel_position_col = int(np.round(pixel_positions[0]))
    pixel_position_row = int(np.round(pixel_positions[1]))

    # 根据所在滤光片或边界情况，计算提取目标区域的半边长（单位：像素）
    # 若 RA,Dec 对应位置接近图像边缘，则做最小防护
    if pixel_position_col < 30 or pixel_position_row < 30:
        half_side_length = min([pixel_position_col, pixel_position_row])
    else:
        if band in ['u', 'v', 'g', 'r']:
            # 滤光片为 u,v,g,r 时，像素分辨率为 0.429 "/pix
            half_side_length = int(target_region_half_size / 0.429)
        elif band in ['i', 'z', 'y']:
            # 滤光片为 i,z,y 时，像素分辨率为 0.286 "/pix
            half_side_length = int(target_region_half_size / 0.286)

    # 从原图像中截取目标区域（2*half_side_length 大小的方形区域）
    target_region = np.array(
        img_data_non_scaled[
            int(pixel_position_row) - half_side_length : int(pixel_position_row) + half_side_length,
            int(pixel_position_col) - half_side_length : int(pixel_position_col) + half_side_length
        ]
    )

    # 为防止字节序问题（如 big-endian / little-endian），转换为原生字节序
    target_region = target_region.byteswap().newbyteorder()

    # 使用 SEP 库对该目标区域进行背景估计，得到背景 RMS
    target_region_bkg_map = sep.Background(target_region, bw=64, bh=64, fw=3, fh=3)
    target_region_bkg_rms_global = target_region_bkg_map.globalrms

    # 再次读取图像头部信息，提取 FWHM（假定已存储在主头里）
    with fits.open(image_path) as hdul:
        header_ = hdul[0].header
    if "arcsec" in header_.comments["FWHM"]:
        median_fwhm = header_["FWHM"] / PIXSCAL  # 假设此关键字中存放了测量得到的 FWHM
        logger.info(f'''header_["FWHM"]:{header_["FWHM"]},Unit:arcsec, PIXSCAL:{PIXSCAL}，median_fwhm:{median_fwhm} ''')
    else:
        median_fwhm = header_["FWHM"]
        logger.info(f'''header_["FWHM"]:{header_["FWHM"]},Unit:Pixel, PIXSCAL:{PIXSCAL}，median_fwhm:{median_fwhm} ''')
    # 根据背景 RMS、曝光时间、FWHM 和零点计算 3σ 和 5σ 的极限星等
    # 公式中：2.266 * (median_fwhm^2) 是用于估计光斑面积的一个常数
    limiting_mag_3sigma = -2.5 * np.log10(
        3 * (target_region_bkg_rms_global / exp_time) * np.sqrt(2.266 * (median_fwhm ** 2))
    ) + mag_zp
    limiting_mag_5sigma = -2.5 * np.log10(
        5 * (target_region_bkg_rms_global / exp_time) * np.sqrt(2.266 * (median_fwhm ** 2))
    ) + mag_zp

    # 绘图可视化目标区域，并在标题中显示背景 RMS、FWHM 和 3σ/5σ 极限星等
    plt.figure()
    plt.title(
        "%s_$%s^{\\prime\\prime}\\times%s^{\\prime\\prime}$\n"
        "Background$_\\mathregular{rms}$:%.3f, FWHM:%.1f pixels\n"
        "Limiting magnitude$_{3\\sigma/5\\sigma}$: %.3f/%.3f mag"
        % (
            image_path.split('/')[-1].split('_sciimg')[0].split('.fits')[0],
            target_region_half_size * 2,
            target_region_half_size * 2,
            target_region_bkg_rms_global,
            median_fwhm,
            limiting_mag_3sigma,
            limiting_mag_5sigma
        )
    )
    plt.imshow(
        target_region,
        interpolation='nearest',
        cmap='gray',
        origin='lower',
        vmin=np.nanmin(target_region),
        vmax=np.nanmax(target_region)
    )
    plt.tick_params(labelbottom=False, labelleft=False)  # 隐藏坐标刻度标签

    # 将可视化的 PNG 图片保存
    plt.savefig(
        output_path.replace('.fits', f'_{filter_}_LimitingMag.png'),
        dpi=300, bbox_inches='tight'
    )
    logger.info(output_path.replace('.fits', f'_{filter_}_LimitingMag.png'))
    plt.close()

    # 将计算结果整理成列表返回
    limiting_mag_info = [
        image_path.split('/')[-1],  # 图片文件名
        band,                       # 滤光片
        exp_time,                   # 曝光时间
        median_fwhm,                # FWHM
        target_region_bkg_rms_global,  # 背景 RMS
        limiting_mag_3sigma,           # 3σ 极限星等
        limiting_mag_5sigma,           # 5σ 极限星等
        mag_zp                         # 零点
    ]
    return limiting_mag_info



def calculate_mid_exposure_time(date_value, exposure_value):
    exposure_start_time = datetime.strptime(date_value, '%Y-%m-%dT%H:%M:%S')
    half_exposure_time = timedelta(seconds=exposure_value / 2)
    mid_exposure_time = exposure_start_time + half_exposure_time
    return mid_exposure_time



def bkg_uniform_e2v(image_matrix):
    scaleList = {}

    # 0) linear fitter
    fitter            = fitting.LinearLSQFitter()
    outlier_fitter    = fitting.FittingWithOutlierRemoval(fitter, sigma_clip, niter=3, sigma=3.0)
    line_init         = models.Linear1D()

    # 1) define the pixel boundary for each gate for correction
    pix_bound = {}
    for i in range(16):
        igate = f"gate{i+1}"
        ii = i//8
        iy0, iy1 = 0+4616*ii, 4616*(ii+1)
        if i<=7:
            ix0, ix1 = 0+1152*i, 1152*(i+1)
        else:
            jj = 16-1-i
            ix0, ix1 = 0+1152*jj, 1152*(jj+1)
        
        pix_bound[igate] = [ix0, ix1, iy0, iy1]
        scaleList[igate] = [ix0, ix1, iy0, iy1, 1.0]

    pix_bound_top    = [   0, 9216, 4616, 4621]
    pix_bound_bottom = [   0, 9216, 4611, 4616]

    # 2) extrac the pixels in each boundary
    bins = np.linspace(0,4616,21, dtype=int)
    nbin = len(bins) - 1
    scale_bottom = [1.0]*8
    for i in reversed(range(1,8)):
        igate = f"gate{i+1}"
        jgate = f"gate{i}"
        ix0, ix1, iy0, iy1 = pix_bound[igate]
        jx0, jx1, jy0, jy1 = pix_bound[jgate]
        ipixs = image_matrix[iy0:iy1, ix0:ix0+5]
        jpixs = image_matrix[jy0:jy1, jx1-5:jx1]

        gi, gj = [], []
        for k in range(nbin):
            kbin0, kbin1 = bins[k], bins[k+1]
            kmean1, kmed1, kstd1 = sigma_clipped_stats(ipixs[kbin0:kbin1, :], stdfunc='mad_std')
            kmean2, kmed2, kstd2 = sigma_clipped_stats(jpixs[kbin0:kbin1, :], stdfunc='mad_std')
            gi += [kmed1]
            gj += [kmed2]
        cij = np.corrcoef(gi, gj)[0,1]
        if abs(cij)>0.1:
            fitted_line, mask = outlier_fitter(line_init, np.array(gi), np.array(gj))
            #print(fitted_line.slope.value, fitted_line.intercept.value)
            slope, interp = fitted_line.slope.value, fitted_line.intercept.value
            image_matrix[iy0:iy1, ix0:]    = slope*image_matrix[iy0:iy1, ix0:] + interp
            
            scale_bottom[i:] = list(np.array(scale_bottom[i:]) * slope)
        #plt.plot(gi, gj, "o", color="black")
        #plt.plot(gi, fitted_line(gi), "--", color="black")
        #plt.show()
    
    scale_top = [1.0]*8
    for i in range(9, 16):
        igate = f"gate{i+1}"
        jgate = f"gate{i}"
        ix0, ix1, iy0, iy1 = pix_bound[igate]
        jx0, jx1, jy0, jy1 = pix_bound[jgate]
        ipixs = image_matrix[iy0:iy1, ix1-5:ix1]
        jpixs = image_matrix[jy0:jy1, jx0:jx0+5]

        gi, gj = [], []
        for k in range(nbin):
            kbin0, kbin1 = bins[k], bins[k+1]
            kmean1, kmed1, kstd1 = sigma_clipped_stats(ipixs[kbin0:kbin1, :], stdfunc='mad_std')
            kmean2, kmed2, kstd2 = sigma_clipped_stats(jpixs[kbin0:kbin1, :], stdfunc='mad_std')
            gi += [kmed1]
            gj += [kmed2]
        cij = np.corrcoef(gi, gj)[0,1]
        if abs(cij)>0.1:
            fitted_line, mask = outlier_fitter(line_init, np.array(gj), np.array(gi))
            #print(fitted_line.slope.value, fitted_line.intercept.value)
            slope, interp = fitted_line.slope.value, fitted_line.intercept.value
            image_matrix[iy0:iy1, jx0:]    = slope*image_matrix[iy0:iy1, jx0:] + interp
            
            scale_top[:i-8] = list(np.array(scale_top[:i-8]) * slope)
        #plt.plot(gj, gi, "o", color="black")
        #plt.plot(gj, fitted_line(gj), "--", color="black")
        #plt.show()

    # 4) linear fitting to find the relationship between the top and bottom gates
    pixs_bottom = image_matrix[pix_bound_bottom[2]:pix_bound_bottom[3], pix_bound_bottom[0]:pix_bound_bottom[1]]
    pixs_top    = image_matrix[pix_bound_top[2]:pix_bound_top[3], pix_bound_top[0]:pix_bound_top[1]]

    bins = np.linspace(0,9216,25, dtype=int)
    nbin = len(bins) - 1
    mb, mt = [], []
    for i in range(nbin):
        ibin0, ibin1 = bins[i], bins[i+1]
        imean1, imed1, istd1 = sigma_clipped_stats(pixs_bottom[:, ibin0:ibin1], stdfunc='mad_std')
        imean2, imed2, istd2 = sigma_clipped_stats(pixs_top[:, ibin0:ibin1], stdfunc='mad_std')
        mb += [imed1]
        mt += [imed2]

    cbt = np.corrcoef(mb, mt)[0,1]
    # print(cbt)
    if abs(cbt)>0.1:
        fitted_line0, mask0 = outlier_fitter(line_init, np.array(mb), np.array(mt))
        slope0, interp0 = fitted_line0.slope.value, fitted_line0.intercept.value
        image_matrix[0:4616, :] = slope0*image_matrix[0:4616,:] + interp0
        
        scale_bottom = list(np.array(scale_bottom)*slope0)

    #plt.plot(mb, mt, "o", color="black")
    #plt.plot(mb, fitted_line0(mb), "--", color="black")
    #plt.plot()
    #plt.show()

    for i in range(16):
        igate = f"gate{i+1}"
        if i<=7:
            scaleList[igate][-1] = scale_bottom[i]
        else:
            scaleList[igate][-1] = scale_top[i-8]
    return image_matrix, scaleList

def bkg_uniform_andor(image_matrix):
    scaleList = {}

    # 0) linear fitter
    fitter            = fitting.LinearLSQFitter()
    outlier_fitter    = fitting.FittingWithOutlierRemoval(fitter, sigma_clip, niter=3, sigma=3.0)
    line_init         = models.Linear1D()

    # 1) define the pixel boundary for each gate for correction
    pix_bound_gate1  = [3067, 3072,    0, 3080]
    pix_bound_gate2  = [3072, 3077,    0, 3080]
    pix_bound_gate3  = [3072, 3077, 3080, 6160]
    pix_bound_gate4  = [3067, 3072, 3080, 6160]
    pix_bound_top    = [   0, 6144, 3080, 3085]
    pix_bound_bottom = [   0, 6144, 3075, 3080]
 
    # 2) extrac the pixels in each boundary
    pixs_gate1 = image_matrix[pix_bound_gate1[2]:pix_bound_gate1[3], pix_bound_gate1[0]:pix_bound_gate1[1]]
    pixs_gate2 = image_matrix[pix_bound_gate2[2]:pix_bound_gate2[3], pix_bound_gate2[0]:pix_bound_gate2[1]]
    
    pixs_gate3 = image_matrix[pix_bound_gate3[2]:pix_bound_gate3[3], pix_bound_gate3[0]:pix_bound_gate3[1]]
    pixs_gate4 = image_matrix[pix_bound_gate4[2]:pix_bound_gate4[3], pix_bound_gate4[0]:pix_bound_gate4[1]]

    # 3) linear fitting to find the relationship between the gates
    bins = np.linspace(0,3080,21, dtype=int)
    nbin = len(bins) - 1
    g1, g2, g3, g4 = [], [], [], []
    for i in range(nbin):
        ibin0, ibin1 = bins[i], bins[i+1]
        imean1, imed1, istd1 = sigma_clipped_stats(pixs_gate1[ibin0:ibin1, :], stdfunc='mad_std')
        imean2, imed2, istd2 = sigma_clipped_stats(pixs_gate2[ibin0:ibin1, :], stdfunc='mad_std')
        imean3, imed3, istd3 = sigma_clipped_stats(pixs_gate3[ibin0:ibin1, :], stdfunc='mad_std')
        imean4, imed4, istd4 = sigma_clipped_stats(pixs_gate4[ibin0:ibin1, :], stdfunc='mad_std')

        g1 += [imed1]
        g2 += [imed2]
        g3 += [imed3]
        g4 += [imed4]
    
    # calculate the correlation
    c12 = np.corrcoef(g1, g2)[0,1]
    c34 = np.corrcoef(g3, g4)[0,1]
    #print(c12, c34)
    slope1, interp1 = 1.0, 0.0
    slope2, interp2 = 1.0, 0.0
    if abs(c12)>0.1:
        fitted_line1, mask1 = outlier_fitter(line_init, np.array(g2), np.array(g1))
        #print(fitted_line1.slope.value, fitted_line1.intercept.value)
        slope1, interp1 = fitted_line1.slope.value, fitted_line1.intercept.value
        image_matrix[0:3080,3072:6144]    = slope1*image_matrix[0:3080,3072:6144] + interp1
    if abs(c34)>0.1:
        fitted_line2, mask2 = outlier_fitter(line_init, np.array(g3), np.array(g4))
        #print(fitted_line2.slope.value, fitted_line2.intercept.value)
        slope2, interp2 = fitted_line2.slope.value, fitted_line2.intercept.value
        image_matrix[3080:6160,3072:6144] = slope2*image_matrix[3080:6160,3072:6144] + interp2

    #plt.plot(g2, g1, "o", color="black")
    #plt.plot(g3, g4, "s", color="red")
    #plt.plot(g2, fitted_line1(g2), "--", color="black")
    #plt.plot(g3, fitted_line2(g3), "--", color="red")
    #plt.plot()
    #plt.show()

    # 4) linear fitting to find the relationship between the top and bottom gates
    pixs_bottom = image_matrix[pix_bound_bottom[2]:pix_bound_bottom[3], pix_bound_bottom[0]:pix_bound_bottom[1]]
    pixs_top    = image_matrix[pix_bound_top[2]:pix_bound_top[3], pix_bound_top[0]:pix_bound_top[1]]
    
    bins = np.linspace(0,6144,33, dtype=int)
    nbin = len(bins) - 1
    mb, mt = [], []
    for i in range(nbin):
        ibin0, ibin1 = bins[i], bins[i+1]
        imean1, imed1, istd1 = sigma_clipped_stats(pixs_bottom[:, ibin0:ibin1], stdfunc='mad_std')
        imean2, imed2, istd2 = sigma_clipped_stats(pixs_top[:, ibin0:ibin1], stdfunc='mad_std')

        mb += [imed1]
        mt += [imed2]
    
    cbt = np.corrcoef(mb, mt)[0,1]
    slope0, interp0 = 1.0, 0.0
    # print(cbt)
    if abs(cbt)>0.1:
        fitted_line0, mask0 = outlier_fitter(line_init, np.array(mb), np.array(mt))
        slope0, interp0 = fitted_line0.slope.value, fitted_line0.intercept.value
        image_matrix[0:3080, :] = slope0*image_matrix[0:3080,:] + interp0

    #plt.plot(mb, mt, "o", color="black")
    #plt.plot(mb, fitted_line0(mb), "--", color="black")
    #plt.plot()
    #plt.show()
    
    scaleList["gate1"] = [   0, 3072,    0, 3080, slope0]
    scaleList["gate2"] = [3072, 6144,    0, 3080, slope0 * slope1]
    scaleList["gate3"] = [3072, 6144, 3080, 6160, slope2]
    scaleList["gate4"] = [   0, 3072, 3080, 6160, 1.0]
    return image_matrix, scaleList

def crossmatch(ra1, dec1, ra2, dec2, aperture=1.0):
    """
    Match two sets of on-sky coordinates to each other.
    I.e., find nearest neighbor of one that's in the other.
    """
    """
    Finds matches in one catalog to another.
    
    Parameters
    ra1 : array-like
          Right Ascension in degrees of the first catalog
    dec1 : array-like
          Declination in degrees of the first catalog (shape of array must match `ra1`)
    ra2 : array-like
          Right Ascension in degrees of the second catalog
    dec2 : array-like
          Declination in degrees of the second catalog (shape of array must match `ra2`)
    aperture : cross-matching aperture, float, default 1.0"
                How close (in arcseconds) a match has to be to count as a match.
    Returns
    -------
    idx1 : int array
           Indecies into the first catalog of the matches. Will never be
           larger than `ra1`/`dec1`.
    idx2 : int array
           Indecies into the second catalog of the matches. Will never be
           larger than `ra1`/`dec1`.
    """
    ra1 = np.array(ra1, copy=False)
    dec1 = np.array(dec1, copy=False)
    ra2 = np.array(ra2, copy=False)
    dec2 = np.array(dec2, copy=False)

    if ra1.shape != dec1.shape:
        raise ValueError('!! ra1 and dec1 do not match!')
    if ra2.shape != dec2.shape:
        raise ValueError('!! ra2 and dec2 do not match!')

    nobj1, nobj2 = len(ra1), len(ra2)
    if nobj1 > nobj2:
        ra1, ra2 = ra2, ra1
        dec1, dec2 = dec2, dec1

    x1, y1, z1 = _spherical_to_cartesian(ra1.ravel(), dec1.ravel())
    coords1 = np.empty((x1.size, 3))
    coords1[:, 0], coords1[:, 1], coords1[:, 2] = x1, y1, z1

    x2, y2, z2 = _spherical_to_cartesian(ra2.ravel(), dec2.ravel())
    coords2 = np.empty((x2.size, 3))
    coords2[:, 0], coords2[:, 1], coords2[:, 2] = x2, y2, z2

    kdt = ckdt(coords2)
    idxs2 = kdt.query(coords1)[1]

    ds = _great_circle_distance(ra1, dec1, ra2[idxs2], dec2[idxs2]) # in arcsecond
    idxs1 = np.arange(ra1.size)

    msk = ds < aperture
    idxs1, idxs2 = idxs1[msk], idxs2[msk]
    ds = ds[msk]
    # Usually, there is duplicate ID in idxs2, here we only keep the one with smaller distance
    dupid = [xx for xx, yy in Counter(idxs2).items() if yy > 1]
    badid = np.array([])
    if len(dupid) > 0:
        for k in dupid:
            kid = np.where(idxs2==k)[0]
            nkid = np.delete(kid,ds[kid].argmin())
            badid = np.append(badid,nkid)
    
    if len(badid)>0: idxs1, idxs2 = np.delete(idxs1,badid), np.delete(idxs2,badid)

    if nobj1 > nobj2:
        newid = np.argsort(idxs2)
        idxs1, idxs2 = idxs2[newid], idxs1[newid]

    return idxs1, idxs2

def _spherical_to_cartesian(ra, dec):
    """
    (Private internal function)
    Inputs in degrees. Outputs x,y,z
    """
    rar = np.radians(ra)
    decr = np.radians(dec)

    x = np.cos(rar) * np.cos(decr)
    y = np.sin(rar) * np.cos(decr)
    z = np.sin(decr)

    return x, y, z

def _great_circle_distance(ra1, dec1, ra2, dec2):
    """
    (Private internal function)
    Returns great ciircle distance. Inputs in degrees.

    Uses vicenty distance formula - a bit slower than others, but
    numerically stable.
    """
    lambs = np.radians(ra1)
    phis = np.radians(dec1)
    lambf = np.radians(ra2)
    phif = np.radians(dec2)

    dlamb = lambf - lambs

    numera = np.cos(phif) * np.sin(dlamb)
    numerb = np.cos(phis)*np.sin(phif) - np.sin(phis)*np.cos(phif)*np.cos(dlamb)
    numer = np.hypot(numera, numerb)
    denom = np.sin(phis)*np.sin(phif) + np.cos(phis)*np.cos(phif)*np.cos(dlamb)
    return np.degrees(np.arctan2(numer, denom))*3600.0 # convert to arcsecond



def wds9reg(x,y,flag=None,radius=5.0,unit="arcsec",color="green",outfile="out.reg"):
    """
    Write ds9 region file.

    Parameters:
    coordinate: 2D array
       coordinate to be written.  It could be image coordinates or RA/Dec.  
       Former, set unit="pixel"
       Later, set unit="arcsec"
    radius: float
       in unit of pixels or arcsec (when 'unit' keyword is set)
    unit: string
        pixel: write region file in unit of pixels
        arcsec (default): write region file in unit of RA and Dec
    color: string
       to specify which color to use.  Default is green
    outfile: string
       name of output region file
    
    Return:
       "outfile": can be read by ds9

    Example:
        pos = [100, 200]
        wds9reg(pos,outfile="pos.reg")
    """
    if not unit in ["arcsec","pixel"]:
        raise ValueError("!! Please set 'unit' as 'arcsec' or 'pixel'")

    fileobj = open(outfile, "w")
    note0 = "# Region file for DS9\n"
    global_pro1 = "global color=%s font='helvetica 10 normal' "%color
    global_pro2 = "select=1 edit=1 move=1 delete=1 include=1 fixed=0 source\n"
    fileobj.write(note0)
    fileobj.write(global_pro1+global_pro2)

    if unit == "arcsec":
        fileobj.write("fk5\n")
        fmt = 'circle(%10.6f,%10.6f,%5.2f")\n'
        if flag is not None: fmt='circle(%10.6f,%10.6f,%5.2f") # text={%d}\n'
    if unit == "pixel":
        fileobj.write("image\n")
        fmt = 'circle(%10.6f,%10.6f,%5.2f)\n'
        if flag is not None: fmt='circle(%10.6f,%10.6f,%5.2f) # text={%d}\n'

    for i in range(len(x)):
        if flag is not None:
            ds9row = fmt%(x[i],y[i],radius,flag[i])
        else:
            ds9row = fmt%(x[i], y[i], radius)
        fileobj.write(ds9row)
    fileobj.close()

    return

def deg2str(ra, dec):
    """
    convert ra&dec in (degree, degree) to (hhmmss, ddmmss),

    Parameters:
    ra, dec: 
       if (ra, dec) is in (deg, deg), float
    
    Return:
       ra  = {0:HH:MM:SS.SS, 
              1:HHMMSS.SS}
       dec = {0:(+/-)DD:MM:SS.SS,
              1:(P/N)DD:MM:SS.SS,
              2:(+/-)DDMMSS.SS,
              3:(P/N)DDMMSS.SS}
    Example:
    deg2str(1.0, 1.0, output_fmt=[0,1])
    """
    rah = ra/15.0
    ram = (rah - int(rah))*60.0
    ras = (ram - int(ram))*60.0

    ra_list = [f"0{int(rah)}"[-2:], f"0{int(ram)}"[-2:], f"0{float(ras):.2f}"[-5:]]
    sra = {0: ":".join(ra_list),
           1: "".join(ra_list)}

    decabs = abs(dec)
    decd   = int(decabs)
    decm   = (decabs - decd)*60.0
    decs   = (decm - int(decm))*60.0
    dec_list = [f"0{int(decd)}"[-2:], f"0{int(decm)}"[-2:], f"0{float(decs):.2f}"[-5:]]
    if dec < 0.0:
        sdec = {0: "-" + ":".join(dec_list),
                1: "N" + ":".join(dec_list),
                2: "-" + "".join(dec_list),
                3: "N" + "".join(dec_list)}
    else:
        sdec = {0: "+" + ":".join(dec_list),
                1: "P" + ":".join(dec_list),
                2: "+" + "".join(dec_list),
                3: "P" + "".join(dec_list)}

    return sra, sdec

def str2deg(ra, dec):
    """
    convert ra&dec in (hhmmss, ddmmss) to (degree, degree)

    Parameters:
    ra, dec: string in (hms, dms)

    Example:
    str2deg("00:00:00.00", "+00:00:00.00")
    str2deg("00:00:00.00", "N00:00:00.00")
    """
    decSign = dec[0]
    
    sra = np.array(ra.split(":"), dtype=float)
    sra = ((sra[-1]/60.0+sra[1])/60.0 + sra[0])*15.0
    
    sdec = np.array(dec[1:].split(":"), dtype=float)
    if decSign in ["-", "N"]:
        sdec = -((sdec[-1]/60.0+sdec[1])/60.0 + abs(sdec[0]))
    elif decSign in ["+", "P"]:
        sdec = (sdec[-1]/60.0+sdec[1])/60.0 + sdec[0]
    else:
        raise ValueError("!!! Give a right dec value")
    
    return sra, sdec


import subprocess

def run_command(cmd):
    """
    运行系统命令并返回结果。

    Args:
    - cmd (str): 要执行的系统命令。

    Returns:
    - dict: 包含命令执行结果的信息，包括:
        - success (bool): 是否执行成功。
        - return_code (int): 命令的退出状态码。
        - stdout (str): 标准输出信息。
        - stderr (str): 标准错误信息。
    """
    try:
        # 执行命令并捕获输出和错误
        result = subprocess.run(
            cmd,
            shell=True,
            capture_output=True,
            text=True,
            check=False
        )
        
        # 返回结果字典
        return {
            "success": result.returncode == 0,  # 判断是否执行成功
            "return_code": result.returncode,   # 返回状态码
            "stdout": result.stdout.strip(),    # 标准输出
            "stderr": result.stderr.strip()     # 标准错误
        }
        
    except Exception as e:
        # 捕获可能的异常并返回失败状态
        return {
            "success": False,
            "return_code": -1,
            "stdout": "",
            "stderr": str(e)
        }